Affine Image Registration Transformation Estimation 
Using a Real Coded Genetic Algorithm with SBX 
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This paper describes the application of a real coded genetic algorithm (GA) to align two 
or more 2-D images by means of image registration. The proposed search strategy is a trans- 
formation parameters-based approach involving the affine transform. The real coded GA uses 
Simulated Binary Crossover (SBX), a parent-centric recombination operator that has shown to 
deliver a good performance in many optimization problems in the continuous domain. In addi- 
tion, we propose a new technique for matching points between a warped and static images by 
using a randomized ordering when visiting the points during the matching procedure. This new 
technique makes the evaluation of the objective function somewhat noisy, but GAs and other 
population-based search algorithms have been shown to cope well with noisy fitness evaluations. 
The results obtained are competitive to those obtained by state-of-the-art classical methods in 
image registration, confirming the usefulness of the proposed noisy objective function and the 
suitability of SBX as a recombination operator for this type of problem. 

Image Registration (IR) is the process of aligning two or more images of the same scene taken at 
different times, from different directions, and/or by different sensors, by finding the best mapping 
function between them [2j [32]. Most research in this area is based on classical algorithms and 
methods, but during the past two decades or so, there has been a growing interest in the application 
of Evolutionary Computation (EC) and other Metaheuristic (MH) methods to solve the problem. 

This paper addresses the design of the mapping function for 2-D image registration using the 
affine transform. In this problem, the goal is to find the best mapping function, also called trans- 
form, that warps a Deformed image (D) in the direction of a Static image (S), based in the images' 
features (e.g. point positions). The objective of the registration procedure is to find the optimum, 
or an acceptable sub-optimal, mapping function. 

To solve this problem we use an elitist real coded genetic algorithm using Simulated Binary 
Crossover and Gaussian mutation. The results obtained are superior to those obtained by a previous 
real coded genetic algorithm on the same synthetic images, and are also very competitive with 
classical state-of-the-art image registration algorithms. 
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The paper is organized as follows. Section [T] describes the IR problem and the two broad 
categories of methods commonly used. It also describes the affine transformation, a commonly 
used linear transformation whose parameters correspond to the decision variables of the optimiza- 
tion problem that is formulated. Section [2] presents a brief review of related work that has been 
proposed to address the IR problem, both with classical as well as with EC and MH methods. Sec- 
tion [3] presents a real coded GA formulation for the IR problem. Section [4] presents and discusses 
experimental results. Finally, Section [5] concludes the paper. 

1 Image Registration 

Image registration has been applied in a large number of research areas, including medical image 
analysis, computer vision and pattern recognition [32]. In all IR problems there are at least two 
images, a Deformed {D) and a Static (S) image, that represent the same object, or scene. IR aims 
to find the best mapping function T to warp D towards S, as shown below, 

W^T(D) . (1) 

W is a warped image which should be as closely shaped to S as possible. IR typically has the four 
following steps [32] : 

1. Feature detection; 

2. Feature matching; 

3. Mapping function design; 

4. Image transformation and re-sampling. 

In order to find the correct transformation, image features such as closed-boundary regions, edges, 
line intersections, corners, and so on, should be extracted (feature detection). These features can 
be used as control points. Correspondences have to be set between the extracted features of S and 
D (feature matching). Then, the type of transforming model has to be defined and its parameters 
estimated (mapping function). Finally, the D image is transformed by means of the mapping 
function (image transformation). 

IR can be seen as a function approximation method [24], and it is an NP-Complete problem [20] . 
The most important aspect of IR is the discovery of the unknown parametric transformation that 
relates the two images. Two different approaches can be found in the literature: 

• Matching-based approaches 

• Transformation parameters-based approaches 

Matching-based approaches conduct a search within the space of possible feature correspon- 
dences (typically point matching) between the two images. Thereafter, a registration transforma- 
tion is obtained based on the correspondence found. In contrast, transformation parameters-based 
approaches perform a direct search in the space of the parameters of the transformation. 

IR methods can also be classified according to the type of models that they allow to transform 
the D image into the S image (sometimes also referred in the literature as the scene and model 
images). Two major types of models are used: linear and non-linear transformations. Linear 
transformations preserve the operations of vector addition and scalar multiplication. The same 
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does not hold for non-linear (or elastic) transformations, which allow local deformations of the 
image. Although, with a few modifications, the proposed approach could be adapted to work with 
polynomial transforms (e.g. quadratic and cubic), in this paper we only deal with linear transfor- 
mation models. Several such models exist and one of the most general is the affine transformation 
which is discussed next. 

1.1 Affine Transform 

The affine transform is a linear transformation that includes the following elementary transforma- 
tions: translation, rotation, scaling, stretching, and shearing |19j . These elementary transforma- 
tions are illustrated in Figure CD 




Translation 



Rotation 



Scaling 



Stretch 



Shearing 



Figure 1: Elementary geometric transforms for a planar surface element used in the affine transform: 
translation, rotation, scaling, stretching, and shearing. 

A geometric operation transforms a given image D into a new image W by modifying the 
coordinates of the image points, as follows: 



D(x,y)^W(x',y') 



(2) 



The original values of image D which are located at (x,y) warp to the new position (x',y') in the 
new image W. To model this process, we first need a mapping function T that is a continuous 
coordinate transform. An affine transformation function works in the 2-D space, thus, the search 
space is: 

T : M 2 => M 2 . (3) 



The mapping function can be redefined as: 



W = T(D) 
T : M 2 =>• M 2 



(4) 



The warped image W(x',y') is often specified as the following two separated functions for the x 
and y components: 



x' = T x (x,y) 

y' = T y( x ,y) ■ 



(5) 
(6) 
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An affine transformation can be expressed by vector addition and matrix multiplication as shown 
in Equation [7J 
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where S is the scaling parameter. By multiplying S with the rotation matrix, Equation [7] can be 
written as: 

(8) 

Finally by using homogeneous coordinates, the affine transformation can be rewritten as Equa- 
tion M 
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The affine transform has six parameters: 9o, 0i, 9%, 9s, 0^, and 9g. 62 and 6 '5 specify the translation 
and 9q, 9i, 63, and 9^ aggregate rotation, scaling, stretching, and shearing. 



2 Related Work 

There is variety of techniques for solving IR problems. This section presents a brief review of 
some of the most important ones, including both classical as well as evolutionary computation and 
metaheuristic based approaches. 



2.1 Classical Methods 

Most of the classical approaches for IR can be found in [321 l2l 15]. Two state-of-the-art approaches 
from this category of methods are Robust Point Matching (TPS-RPM) [4j, and Shape Context 

(sc) p. 

TPS-RPM is a method for matching two point-sets in a Deterministic Annealing (DA) setting. 
It uses a fuzzy-like matrix instead of a binary permutation matrix to find the matching between two 
sets of points. In TPS-RPM, both the point correspondences and the transformations are computed 
interchangeably. Therefore, RPM can be viewed as a general framework for point matching and can 
accept different transformation models [3] like affine, and even more complicated models like Thin 
Plate Splines (TPS) [12]. This method is a kind of a hybrid in the sense that it can be considered 
both a matching-based and a transformation-based approach for IR. 

Shape Context (SC) is a matching-based approach that is usually used to estimate the transfor- 
mation between two images, by finding matches between samples from the edges of the objects in 
the images. It basically consists of analyzing the spacial relationship between points. It uses four 
main parameters. The first defines the number of radial bins for the creation of the histograms, the 
second is the number of theta bins that defines how many slices should the histograms be divided 
into, and the third and fourth parameters, the minimum and maximum width of the bins. For 
more information on these and other classical IR methods, the reader is directed to [2~1 \5\ [32]. 



2.2 EC and MH methods 

Evolutionary Algorithms (EAs) and other metaheuristics (MHs) have been applied to solve IR 
problems. EAs and MHs are stochastic optimization methods whose goal is to find a solution or a 
set of solutions that perform(s) best with respect to a certain objective(s). During the last decades 
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these algorithms have been successful in solving a variety of search and optimization problems, 
and the domain of image registration has been no exception. As opposed to the classical methods, 
which are typically based on gradient-based search, EAs and MHs tend to escape more easily from 
local optima and can be considered, in general, robust methods. 

The first known application of evolutionary computation to image registration is due to Fitz- 
patrick et al. [13] who applied a genetic algorithm (GA) to relate angiographic images. For the 
subsequent 15 years or so, other EC approaches have been proposed by different authors, but most 
of them were based on the canonical GA with proportionate selection and a binary representation 
for solutions. Such a GA has severe limitations when solving optimization problems in the contin- 
uous domain, especially due to the problem of Hamming cliffs originated from the discretization of 
real valued variables into binary coded values, to the fixed precision that depends on the number 
of bits used for each decision variable, and for imposing lower and upper bounds for a variable's 
value. Moreover, it is known for several years that fitness proportionate selection methods have 
several drawbacks when compared to ordinal-based selection methods such as ranking, tournament, 
or truncation selection |16| . Nonetheless, most of the early EC approaches for IR used such kind of 
GA setup [25\ [26l [T4"l [30], [3T] . Another limitation of the early approaches was that they only dealt 
with translation and rotation |13 [ 114 1 Il8 ] 130]. ignoring scaling, stretching, and shearing. 

Most modern EC applications to IR use a direct real coded representation of solutions \21\ 
HHJ CEl El El E9j [24]- Besides EC, other MH approaches have been applied to IR, namely Tabu 
Search [27], Particle Swarm Optimization [28], Iterated Local Search [8], and Scatter Search [71 [23] 
to name a few. A detailed review of these works cannot be made in the paper, but the interested 
reader can consult recent surveys on the topic [21 [22] . 

3 A real coded GA with SBX for Image Registration 

This section introduces a real coded genetic algorithm for the optimization of the parameters of 
an affine transformation for the case of 2-D images. The proposed algorithm is a transformation 
parameters-based approach, since we are performing a direct search for the parameters that define 
the registration transformation. For the sake of simplicity, we assume we have two 2-D synthetic 
point-sets representing features from the two images. In other words, it is assumed that the feature 
detection step of the IR pipeline has been solved beforehand. In the next subsections, we describe 
the representation, the operators, and the objective function that was used in this study. 

3.1 Representation and Operators 

The representation is straightforward. For the 2-D case, the affine transformation is defined by 
six parameters, Qq . . .65, as explained in Section [1.11 A candidate solution for the GA is therefore 
represented by a chromosome vector with six genes, each a real number. 

With respect to the GA variation operators, we use Simulated Binary Crossover (SBX) proposed 
by Deb and Agrawal [10] and Gaussian mutation. The utilization of SBX crossover and Gaussian 
mutation is a natural choice because the problem has a continuous search space. SBX uses a 
probability distribution to create the offspring, and it does so by biasing the offspring to be created 
near the parents. SBX is a parent- centric recombination operator because the offspring it produces 
are located around the parents. This behaviour contrast with mean-centric recombination operators 
whose offspring are located at the center of mass of parents. Examples of the latter category are 
the unimodal normal distribution crossover (UNDX), simplex crossover (SPX), and blend crossover 
(BLX). 
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It has been shown that parent- centric operators have in general a better performance than mean- 
centric operators This has motivated our choice of SBX as a crossover operator. Surprisingly, 
none of the real coded GAs proposed in the literature for addressing the IR problem has used SBX 
or other parent-centric crossover operators. Instead most used recombination operators such as 
uniform [2TJ [3] , arithmetic [18] and blend crossover [6] . 

3.2 Objective function 

In order to guide the search for an appropriate set of parameters for the affine transformation, 
we need to measure the proximity between the static and the warped image (the deformed image 
after the affine transformation is performed). The closer the two images are, the better the affine 
transformation is. Since each image is represented by a point-set, we need a way to find the 
similarity between two point-sets. To do so we first find a correspondence between points in the 
warped and static images. Once the correspondence is obtained, the objective function value is the 
weighted similarity of the two point-sets using the Euclidean distance of the matched points. 

Algorithm [TJ gives details of the steps involved in evaluating a candidate solution. In the 
algorithm, upper case letters denote matrices or vectors, and lower case letters with subscripts 
denote a specific element of the matrix or vector. The next paragraphs describe the major steps of 
Algorithm [TJ 

To compute the objective function value of a candidate solution, we start by warping the de- 
formed image D according to the parameters of the affine transformation specified in the candidate 
solution C, yielding a new point-set W (line 1 of Algorithm [T]) . Then the matching of points in 
W into points in S is modeled using a correspondence binary matrix M(n, k) based on the closest- 
point rule, where n and k correspond to the number of points in the warped and static images, 
respectively. The closest-point is measured using the Euclidean distance between matched points. 
Each point in any set corresponds, at most, to one point in the other set. To find the correspon- 
dence for each point, the closest point in the other set is chosen. If the nearest point has already 
been assigned to another point, the next non-assigned nearest point is chosen. This procedure is 
performed once for finding the correspondence matrix M'(n,k) from the warped set to the static 
set (W — > S), and another time for finding the correspondence matrix M"(n, k) from the static set 
to the warped set {S — > W). This is achieved by lines 3-24 of Algorithm [TJ 

The order in which the correspondence points are found (lines 17 and 21 of Algorithm [TJ) , plays 
a vital role in the resulting correspondence matrix. A match-order vector is proposed to specify 
the order in which the points of a given set are visited when finding its closest-point match from 
the other set (W — > S and S — > W). For different orderings, different correspondences may be 
found. Therefore, the match-order vector is randomly created in each generation. This makes the 
evaluation of a candidate solution a somewhat noisy process. In a given generation, two identical 
solutions obtain the same objective function value. But the same thing is not necessarily true for 
two identical solutions from different generations. Fortunately, GAs are well known for being able 
to handle well noisy fitness evaluations due to processing a population of solutions. Note that we 
could have used a fixed pre-determined ordering for all evaluations but decided not to do so because 
the point matching procedure would be somewhat biased with respect to the used ordering. 

At this point (line 24 of Algorithm [TJ) , matrices M' and M" specify the point-matchings from 
W — > S and S — >• W, respectively. We then obtain matrices Q and Q* . Q is simply the sum of 
M' and M", thus each qij can have a value of 2, 1, or 0, depending on whether point i matches 
point j in both directions, in a single direction, or has no match at all. Matrix Q* is obtained 
from Q by inverting the non-zero elements. Finally, matrix M is calculated from the element-wise 
multiplication of matrices M' and Q*. The objective function is based on the weighed similarity of 



6 



Algorithm 1 Objective function 



Input: S,D,C 

/* S is the static image points */ 
/* D is the deformed image points */ 
/* C is a chromosome */ 
Output: Objective function value 
W T(D,C); 

/* Euclidean distance between the point-sets */ 
for all i G {1, . . . , n} do 
for all j e {l,...,fc} do 

Sij <- |K,Sj||; 
end for 
end for 

/* Initialize correspondence matrices */ 
for all i € {1, ... , n} do 
for all j E {1, . . . , k} do 

m{ • «- 0; 

end for 
end for 

/* Find the closest non-assigned point */ 
/* O is the matched-order vector */ 
for all i G O do 

j <-W2S(A, M'.i); //(W->S) 
M[. <- 1; 
end for 

for all f S O do 

i^S2W(A,M",j); //(S^W) 
M£<-1; 
end for 
Q = M' + M"; 
/* Calculate weights */ 
for all i G {1, . . . , n} do 
for all j G {l,...,k} do 
if ^ then 

else 

9« <- 0; 
end if 
end for 
end for 

/* Weighting matches */ 
for all i G {1, . . . , n) do 
for all j G {1, . . . ,k} do 

m ij <~ m'ij x <?*■; 
end for 
end for 
fitness <— 0; 
for all i G {1, . . . , n) do 
for all j G {1, . . . ,k} do 

fitness <— fitness + (m^ x <5y); 
end for 
end for 

return fitness; 
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two point-sets using the Euclidean distance of the matched points. The points that are connected 
exclusively from one direction (either W — > S or S —¥ W) are penalized, and those that are 
connected in both directions are given half weight in terms of Euclidean distance. In other words, 
if the connection exists in both directions the objective function value decreases. 

This objective function is very similar to the one used by Seixas et al. [24J. The main difference 
is the use of a newly generated matched order vector in each generation, which makes the point- 
matching procedure less dependent on a fixed ordering of visiting the points. 

4 Experimental Results 

This section describes the experimental results from testing the proposed GA formulation. Five 
point-sets available at |http: //noodle. med.yale.edu/~chui/rpm/T PS-RPM. zip| are used. Each 
set is composed at most by 105 points. They include the deformed and static points' locations. 
The deformed points were generated from the static ones by a non-affined (i.e. free-form) transform. 
This means that it will not be possible to obtain a perfect matching of the images by using an affine 
transformation model alone. All point coordinate values have a precision of 16 floating point. 

The GA setup was the same for all data sets. Most parameter settings were tuned beforehand, 
and held fixed for all the experiments. We use tournament selection without replacement of size 5, 
SBX crossover with distribution index 2 [ID], and Gaussian mutation for all the genes. The crossover 
probability was set to 1.0 and each gene undergoes SBX with probability 0.5. For replacement we 
use a replace worst strategy, with the worst half of the individuals of the current population being 
replaced by the best half of the newly generated solutions. This replacement strategy makes the GA 
elitist, never losing the best solution found so far. The GA ran for 500 iterations. The experiments 
were performed with populations of size 30, 60, 120, 240, and 480 individuals, and for each size, 
100 independent runs were executed. 

In accordance to population sizing theory of GAs [17] , larger populations sizes tend to produce 
a better solution quality, but also at the expense of more processing time. Figure [3] shows the 
objective function value of the best individual in the population at every generation, averaged over 
the 100 runs, for the various populations sizes and for the various point-sets. The performance 
behavior is more or less identical for all images. We can observe a substantial progress for the first 
50 generations, still some progress between generations 50-200, and from there on the improvements 
are minor. As expected, larger population sizes gives better solution quality but the improvements 
are negligible for populations sizes larger than 120. 

Figure [2] illustrates the five point-sets before and after warping. It should be noticed that these 
are sets where deformations are non-affine, therefore the resulting warped images will never match 
perfectly. Nevertheless, they present a very good approximation. When compared to the approach 
from |24] it is possible to observe that our results are more precise. 

In order to assess the quality of the proposed approach in what concerns affine deformations, 
affine deformed images were generated from the same five static images, and used in subsequent 
experiments. The results were compared to those produced by well known classical state-of-the- 
art approaches such as SC and TPS-RPM, and are illustrated in Figure [H It can be seen that 
the results produced by the proposed GA are slightly better than those of the SC and are very 
competitive with those obtained by TPS-RPM. 
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Figure 2: Non-affine distorted point-sets (left, blue dots for static image points, red dots for de- 
formed image points) and GA affine image registration (right, red dots are the warped image points) 
results obtained after 500 generations using population size 120. The warped images are zoomed 
for better visualization. Note that even better matching could be obtained with larger population 
sizes, but the improvements are negligible as shown in Figure [3) 
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Figure 3: Best objective function value through generations for various population sizes obtained 
for various point-sets. The results are averaged over 100 independent runs. 



5 Summary and Conclusions 

This paper proposed a real coded GA that is especially suited for doing image registration of affine 
distorted images. As opposed to previous EC approaches for solving IR, our method uses Simulated 
Binary Crossover, a parent-centric recombination operator that has been giving good results on a 
variety of continuous real world optimization problems within a GA framework. The use of a 
randomized ordering when visiting points during the point-matching procedure was also proposed, 
and although this technique yields a noisy fitness function evaluation, the results obtained show 
that the GA is capable of dealing with it quite well. 

The resulting algorithm was applied to 2-D synthetic point-sets, with deformed images points 
obtained from both affine and non-affine transformations. For the case of non-affine distorted points, 
our method produces a more precise registration than previously published results by means of an 
evolutionary algorithm on the same point-sets [24]. For the case of affine distorted points, the 
proposed real coded GA produced better results than SC, and was competitive with TPS-RPM, 
two well known classical state-of-the-art image registration methods. 
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points are almost on top of each other, meaning that the match is almost perfect. For SC the 
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